//
// Author 
// Xianguo Lu
// lu@physi.uni-heidelberg.de 
//

#ifndef HISTBASE_H
#define HISTBASE_H

#ifndef MATHBASE_H
#include "MathBase.h"
#endif

class HistBase
{
 public:
  enum{
    KPNG = 1,
    KPDF = 2,
    KEPS = 4,
    KROOT = 8,
    KC = 16,
    KJPG = 32, 
    KNORX = 0,
    KNORY = 1
  };

  HistBase(){}
  ~HistBase(){  }
  static void IniDraw(TPad *cc, const Double_t rm=0.02, const Double_t lm = 0.125);
  static void PrintOut(const TString pn, const TCanvas *cc=0, const Int_t kopt=1);
  static void PrintOut(Int_t &id, const TString pn, const TCanvas *cc=0, const Int_t kopt=1);

  static TH1D *ProfileX(const TH2 *hh, const Bool_t kbino = kFALSE, const TString tag="", FILE *flog=0, Double_t *rmserr=0x0);
  static TH1D *BinomialDivide(const TH1 *h1, const TH1 *h2, const TString hname);
  static TH1D * ProjectionX(const TH2 * h2d, const Int_t yb0=0, Int_t yb1=10000000, const Double_t fss = -999);

  static TH1D* GetSliceX(const TH2D *h2d, const Double_t xx, const TString hname ="pj");
  static void MultiNor(const TH2 * h2d, const Int_t nh,  TH1D * h1s[], const Int_t inor, Double_t &hmin, Double_t &hmax);
  static TGraphPolar * ToPolar(TCanvas *cc, const TH1 * hh, const Bool_t kdeg2rad, const Double_t phi0);
  static Bool_t ToPDF(TH1 *hist, const Bool_t kpr=1);
  static void ToGraph(const TString fintxt, const Int_t istart, const Int_t ny, TGraphAsymmErrors * gr[], const Double_t norix[][4], const Bool_t kpos, const Bool_t kconvert, const Double_t *ymax=0x0, const Double_t *ymin=0x0, const Bool_t cutyerr=kFALSE, const Bool_t kxerr=kTRUE);
  static void NormalizeSlicesX(TH2D * hh, const Int_t nthres=1);

  static TH1D * DrawFunc(const TH1 * hh, Double_t (*func)(const Double_t *, const Double_t *), const Double_t *par, const TString opt = "hist C", const Int_t color=kRed);
  static TGraph * DrawFunc(const TGraph * gr, Double_t (*func)(const Double_t *, const Double_t *), const Double_t *par, const Int_t color=kRed,  const Int_t npch0=0, const Double_t *patch0=0x0, const Int_t npch1=0, const Double_t *patch1=0x0);

  static void DrawPull(TPad * & ppull, TH1D * & hpull, TF1 * &f0, const TH1 * hh, const TH1 *hfit, const Double_t x0=0.37, const Double_t y0=0.65, const Double_t x1=0.7, const Double_t y1=0.88);
  static void DrawPar(TPad * &ppar,  TPaveText * &pt, const TObject * hh, const TString sfunc, const Int_t npar, const TString spar[], const Double_t par[], const Double_t err[], const Bool_t kpos=kFALSE, const Double_t px=0, const Double_t py=0);

  static TGraphAsymmErrors * Hist2Graph(const TH1 *hist, const Bool_t delog=kFALSE, const TString xtit="");
  static TGraphErrors *ProjectRMS(const TH2D * h2, const Bool_t xlog=kFALSE);

  protected:
  static void SetColor(TPad *cc);
};

void HistBase::IniDraw(TPad *cc, const Double_t rm, const Double_t lm)
{
  gStyle->SetPalette(1);
  gStyle->SetOptTitle(1);
  gStyle->SetOptStat(0);
  gStyle->SetTitleBorderSize(-1);
  //gStyle->SetTitleX(0.2);
  gStyle->SetTitleX(0.1);
  gStyle->SetTitleW(0.9);
  gStyle->SetTitleStyle(0);
  
  TGaxis::SetMaxDigits(3);

  gPad->SetGrid();
  if(lm>0)
    gPad->SetLeftMargin(lm);
  if(rm>0)
    gPad->SetRightMargin(rm);
  gPad->SetTopMargin(0.08);
  gPad->SetTicks(1,1);
  gPad->Range(0,0,1,1);

  SetColor(cc);
}

void HistBase::SetColor(TPad *cc)
{
  gStyle->SetHistFillColor(0);
  //gStyle->SetFillColor(0);//it conflicts with color palette
  gStyle->SetFrameFillColor(0);
  gStyle->SetPadColor(0);
  gStyle->SetCanvasColor(0);
  gStyle->SetStatColor(0);
  gStyle->SetTitleFillColor(0);

  cc->SetFillColor(0);//this is the desired one!!!

  const Int_t nRGBs = 5;
  const Int_t nCont = 100;
  Double_t stops[nRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  Double_t red[nRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  Double_t green[nRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  Double_t blue[nRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  const Int_t cgc = TColor::CreateGradientColorTable(nRGBs, stops, red, green, blue, nCont);      
  Int_t colors[nCont];
  gStyle->SetNumberContours(nCont); 
  for(Int_t ii=0; ii<nCont; ii++) 
    colors[ii]=cgc+ii;
  gStyle->SetPalette(nCont, colors); 
}

void HistBase::PrintOut(Int_t & id, const TString pn, const TCanvas *cc, const Int_t kopt)
{
  PrintOut(Form("ID%03d_%s", id++, pn.Data()), cc, kopt);
}

void HistBase::PrintOut(const TString pn, const TCanvas *cc, const Int_t kopt)
{
  if(!kopt)
    return;

  if(kopt & KPNG)
    cc->Print(pn+".png");
  if(kopt & KPDF)
    cc->Print(pn+".pdf");
  if(kopt & KROOT)
    cc->Print(pn+".root");
  if(kopt & KEPS)
    cc->Print(pn+".eps");
  if(kopt & KC)
    cc->Print(pn+".C");
  if(kopt & KJPG)
    cc->Print(pn+".jpg");
}

TH1D* HistBase::ProfileX(const TH2 *hin, const Bool_t kbino, const TString tag, FILE *flog, Double_t *rmserr)
{
  const Double_t eps = 1e-10;
  const TString hname = hin->GetName();
  const Double_t xmin = hin->GetXaxis()->GetXmin();
  const Double_t xmax = hin->GetXaxis()->GetXmax();
  const Int_t nx      = hin->GetNbinsX();
  TH1D *hout=new TH1D(hname+"_proj","",nx, xmin, xmax);

  for(Int_t ix=1; ix<= nx; ix++){
    TH1D *calcb = new TH1D(Form("calcb%d",ix),"",1, hin->GetYaxis()->GetXmin(), hin->GetYaxis()->GetXmax());

    Double_t ntot = 0;
    for(Int_t iy = 1; iy <= hin->GetNbinsY(); iy++){
      const Double_t tmpy = hin->GetYaxis()->GetBinCenter(iy);
      const Double_t tmpn = hin->GetBinContent(ix, iy);

      ntot += tmpn;
      if(calcb->Fill(tmpy, tmpn)<0){
        printf("HistBase::ProfileX calcb range fail!! %f %f -- %f %f\n", hin->GetYaxis()->GetXmin(), hin->GetYaxis()->GetXmax(), tmpy, tmpn);
        exit(1);
      }
    }
    if(ntot==0){
      delete calcb; 
      continue;
    }

    const Double_t mean = calcb->GetMean();
    Double_t rms = calcb->GetRMS();
    if(kbino)
      rms = rms/sqrt(ntot);
    if(rms==0)
      rms = eps;
    if(rmserr){
      rmserr[ix] = calcb->GetRMSError();
    }

    if(tag!=""){
      char outlog[500];
      sprintf(outlog, "%s %4.1f %6.0f %6.0f %10.3f +/- %10.3f\n", tag.Data(), hin->GetXaxis()->GetBinLowEdge(ix), hin->GetBinContent(ix, 2), hin->Integral(ix, ix, 1, 2), mean, rms);
      printf(outlog);
      cout<<outlog<<endl;
      if(flog) fprintf(flog, outlog);
    }
    hout->SetBinContent(ix, mean);
    //hout->SetBinError(ix, rms/2.);
    hout->SetBinError(ix, rms);//as root does Sqrt(N)

    delete calcb;
  }

  hout->SetYTitle(hin->GetYaxis()->GetTitle());
  hout->SetXTitle(hin->GetXaxis()->GetTitle());
  hout->GetYaxis()->SetTitleOffset(hin->GetYaxis()->GetTitleOffset());
  hout->SetTitle(hin->GetTitle());
  
  return hout;
}


TH1D *HistBase::BinomialDivide(const TH1 *h1, const TH1 *h2, const TString hname)
{
  const Double_t eps = 1e-10;

  const Double_t xmin = h1->GetXaxis()->GetXmin();
  const Double_t xmax = h1->GetXaxis()->GetXmax();
  const Int_t nx = h1->GetNbinsX();
  TH1D *hout = new TH1D(hname, "", nx, xmin, xmax);

  Int_t nent=0;
  for(Int_t ii=1; ii<=nx; ii++){
    const Double_t n1 = h1->GetBinContent(ii);
    const Double_t n2 = h2->GetBinContent(ii);
    if(n1<0 || n2<0){
      printf("error negative in HistBase::BinomialDivide!! %f %f %f\n", h1->GetBinCenter(ii), n1, n2);
      exit(1);
    }
    if(n1>n2){
      printf("warning n1 > n2 in HistBase::BinomialDivide!! %f %f %f\n", h1->GetBinCenter(ii),n1, n2);
      continue;
    }
    if(n2==0)
      continue;
    
    const Double_t yy = n1/n2;

    Double_t err = sqrt(yy*(1-yy)/n2);    
    if(err==0) 
      err = eps;

    hout->SetBinContent(ii, yy);
    hout->SetBinError(ii, err);
    nent++;
    //printf("test %f %f -- %f %f\n",n1, n2, yy, err);
  }

  hout->SetXTitle(h1->GetXaxis()->GetTitle());
  hout->SetEntries(nent);
  return hout;
}

TH1D * HistBase::ProjectionX(const TH2 * h2d, const Int_t yb0, Int_t yb1, const Double_t fss)
{
  const TAxis * ax = h2d->GetXaxis();
  const TAxis * ay = h2d->GetYaxis();
  if(yb1 > ay->GetNbins()+1)
    yb1 = ay->GetNbins()+1;

  const Int_t nx = ax->GetNbins();
  const Double_t x0 = ax->GetXmin();
  const Double_t x1 = ax->GetXmax();
  const TString xtit = ax->GetTitle();
  TH1D * hout = new TH1D(Form("%sprojx", h2d->GetName()), "", nx, x0, x1);
  hout->SetXTitle(xtit);

  Double_t totnpp = 0, totncc = 0;
  for(Int_t ix = 0; ix <= nx+1; ix++){
    Double_t cont = 0;
    Double_t prescl = -999;
    for(Int_t iy = yb0; iy <= yb1; iy++){
      const Double_t ent = h2d->GetBinContent(ix, iy);
      const Double_t err = h2d->GetBinError(ix,iy);
      if(ent==0)
        continue;

      //ent = n/s, err = sqrt(n)/s
      if(prescl==-999){
        prescl = ent/err/err;
        //printf("HistBase::ProjectionX pre-sclae found: %d %f %f\n", ix, ax->GetBinCenter(ix), prescl);
      }
      else if(prescl-ent/err/err >1e-6){
        printf("HistBase::ProjectionX scale not consistent %f %f -- %f %f\n", prescl, ent/err/err, ent, err);
        exit(1);
      }

      cont += ent;
    }
    if(cont==0)
      continue;

    totnpp += prescl;
    totncc += cont*prescl;

    //negative binomial: varicance-t = e(1-p)/p^2; 
    //e = prescl, p=1/cont;
    //err of cont: sqrt(variance-t)/e = sqrt((1-p)/p^2/e)

    const Double_t sigma = fss>0 ? fss*cont : sqrt( (1-1/cont)/pow(1/cont,2)/prescl );

    //printf("testerror %f %f -- %f -- %f\n", prescl, cont, sigma, sqrt(cont));

    hout->SetBinContent(ix, cont);
    hout->SetBinError(ix, sigma);
  }

  printf("HistBase::ProjectionX NPrescale %f NCount %f\n", totnpp, totncc);

  hout->SetMinimum(0);
  return hout;
}

//better use bin center as xx
TH1D* HistBase::GetSliceX(const TH2D *h2d, const Double_t xx, const TString hname)
{
  const Int_t ix = h2d->GetXaxis()->FindBin(xx);

  TAxis * ay = h2d->GetYaxis();

  const Int_t nbiny = ay->GetNbins();
  const Double_t ymin = ay->GetXmin();
  const Double_t ymax = ay->GetXmax();

  const TString h1n=Form("%s%s", h2d->GetName(), hname.Data());
  TH1D * hout=0;
  gDirectory->GetObject(h1n,hout);
  delete hout;
  hout=new TH1D(h1n,"",nbiny, ymin, ymax);

  Double_t tot=0;
  for(Int_t iy= 0; iy<= nbiny+1; iy++){//including underflow and overflow
    const Int_t cy = hout->FindBin(ay->GetBinCenter(iy));
    
    const Double_t ent = h2d->GetBinContent(ix,iy);
    const Double_t err = h2d->GetBinError(ix,iy);

    hout->SetBinContent(cy, ent);
    hout->SetBinError(cy, err);
    tot+=ent;
  }   
  hout->SetEntries(tot);
 
  TString xtit=h2d->GetXaxis()->GetTitle();
  char tmp[500];
  sprintf(tmp," in (%.2f, %.2f) ",h2d->GetXaxis()->GetBinLowEdge(ix),h2d->GetXaxis()->GetBinUpEdge(ix));
  if(xtit.Contains("[")){
    xtit.Insert(xtit.First("["),tmp);
  }
  else
    xtit.Append(tmp);

  hout->SetTitle(xtit);  
  hout->SetXTitle(ay->GetTitle());
  hout->SetMinimum(0);

  return hout;
}

void HistBase::MultiNor(const TH2 * h2d, const Int_t nh,  TH1D * h1s[], const Int_t inor, Double_t &hmin, Double_t &hmax)
{
  const Int_t nbx = h2d->GetNbinsX();
  const Int_t nby = h2d->GetNbinsY();

  for(Int_t ii=0; ii<nh; ii++){
    const TString tmphn=Form("hMultiNor%d",ii);
    gDirectory->GetObject(tmphn,h1s[ii]);
    delete h1s[ii];
    h1s[ii]=new TH1D(tmphn,"", nbx, h2d->GetXaxis()->GetXmin(), h2d->GetXaxis()->GetXmax());
  }

  hmin = 1e10;
  hmax = -1e10;
  for(Int_t ix=h2d->GetXaxis()->GetFirst(); ix<=h2d->GetXaxis()->GetLast(); ix++){
    Double_t ntot = -999;
    if(inor==KNORY){
      ntot = h2d->Integral(ix,ix, 0, nby+1);
      if(ntot<1)
        continue;
    }
    for(Int_t iy=1; iy<=nby; iy++){
      if(inor==KNORX){
        ntot = h2d->Integral(0, nbx+1, iy, iy);
        if(ntot<1)
          continue;
      }
      const Double_t ent = h2d->GetBinContent(ix, iy);
      const Double_t fra= ent/ntot;
      const Double_t err = sqrt(fra*(1-fra)/ntot);

      if(fra>0){
        if(fra>hmax) hmax = fra;
        if(fra<hmin) hmin = fra;
      }

      h1s[iy]->SetBinContent(ix, fra);
      h1s[iy]->SetBinError(ix, err);
      //printf("%d %d %.1e %.1e %.1e %.1e\n", ix, iy, ent, ntot, fra, err);
      h1s[iy]->SetEntries(ntot);
    }
  }
}

TGraphPolar *  HistBase::ToPolar(TCanvas *cc, const TH1 * hh, const Bool_t kdeg2rad, const Double_t phi0) 
{
  const TAxis * ax = hh->GetXaxis();
  const Int_t nx = ax->GetNbins();

  Double_t xx[nx];
  Double_t yy[nx];
  Double_t ey[nx];

  Int_t npoint=0;

  for(Int_t ix=1; ix<=nx; ix++){
    xx[npoint] = ax->GetBinCenter(ix);
    if(kdeg2rad)
      xx[npoint] *= TMath::DegToRad();

    if(xx[npoint]< phi0 || xx[npoint]> phi0+TMath::TwoPi()){
      printf("HistBase::ToPolar %d %f excluded, phi0 %f\n", ix, xx[npoint], phi0);
      continue;
    }

    yy[npoint] = hh->GetBinContent(ix);
    ey[npoint] = hh->GetBinError(ix);
    npoint++;
  }

  TGraphPolar * gr = new TGraphPolar(npoint, xx, yy, 0, ey);
  gr->SetTitle(hh->GetTitle());

  gr->SetMarkerSize(0.8);
  gr->SetMarkerStyle(20);
  gr->SetMarkerColor(kBlue);
  gr->SetLineColor(kRed);
  gr->SetLineWidth(3);

  gr->Draw("AEP");

  // Update, otherwise GetPolargram returns 0; important for GetPolargram()
  cc->Update();

  gr->GetPolargram()->SetToRadian();

  gr->GetPolargram()->SetNdivRadial(505);
  gr->GetPolargram()->SetNdivPolar(106);

  return gr;
}

Bool_t HistBase::ToPDF(TH1 *hist, const Bool_t kpr)
{
  if(hist->GetEntries()==0)
    return kFALSE;

  Int_t const nbin = hist->GetNbinsX();

  Double_t ntot = 0;
  //check if bin content is number of entries
  for(Int_t ib=0; ib<=nbin+1; ib++){
    const Double_t cont = hist->GetBinContent(ib);
    const Double_t err = hist->GetBinError(ib);
    if(cont==0)
      continue;

    if(err==0){
      printf("incorrect err in HistBase::ToPDF!!\n");
      exit(1);
    }

    //cont = N/C, err = sqrt(N)/C
    const Double_t val=pow(cont/err,2);

    const Int_t ival= (Int_t) (val+0.5);
    if(fabs(val-ival)>1e-2){
      printf("incorrect guess in HistBase::ToPDF! ib/nb: %d/%d cont: %e err %e val %e ival %d d %e\n",ib, nbin, cont, err, val, ival,fabs(val-ival));
      exit(1);
    }
    
    hist->SetBinContent(ib, ival);

    ntot += ival;
  }
   
  //Warning in <TH1F::Sumw2>: Sum of squares of weights structure already created
  //hist->Sumw2();

  //not correct!!
  //hist->Scale(1./totn);

  //now normlaize
  //pdf * bin_width = bin_content / totn = eff 
  //pdf = bin_content / totn / bin_width
  //==> eff and pdf independent of rescaling
    
  //d(pdf*bin_width)=sqrt(eff * (1 - eff) / totn )
  //d pdf = sqrt(eff * (1 - eff) / totn ) / bin_width 
  //d pdf depends on rescaling via totn !! So the histogram must NOT have been scaled!!
  //d pdf/pdf = d bin_content / bin_content
  //          = sqrt ( (1-eff)/ bin_content )
  //          = sqrt ( (1/eff-1)/total_N   )
  //assuming eff<0.5
  //d pdf / pdf > 1 / sqrt( total_N ), with total_N = 1e6, dpdf/pdf > 1e-3!!
  //==> For a data set with statistics N, the relative bin content accuracy is larger than 1/sqrt(N), so normally 1e-3 or worse. Value estimation of numerical function just needs a corresponding (or roughtly higher) accuracy.
  //==> For data with very large N and the numerical accuracy is low, the error in evaluating the function should be taken into account in fitting, i.e. added to chi2.
    
  const Double_t bw = hist->GetBinWidth(1);
  if(kpr)
    printf("======> HistBase::ToPDF: totn: %e\tBin Width: %e\n", ntot, bw);
    
  if(ntot==0){
    printf("HistBase::ToPDF ntot==0 ! %f %d\n", hist->GetEntries(), nbin);
    exit(1);
  }
    
  for(Int_t ib=0; ib<=nbin+1; ib++){
    const Double_t cont = hist->GetBinContent(ib);

    //in case of finit number of bins (i.e. eff not always small), Binomial error is more accurate than Poisson error
    const Double_t eff = cont/ntot;
    const Double_t pdf = eff/bw;

    //set error for 0-bin, as if cont = 1
    const Double_t dpdf = cont==0? sqrt(1-1/ntot)/ntot/bw : sqrt(eff*(1-eff)/ntot) / bw;
    hist->SetBinContent(ib, pdf);
    hist->SetBinError(ib, dpdf);
  }
  hist->SetEntries(ntot);

  return kFALSE;
}

//istart from 0
//kconvert = 1, x = log10bg converts to bg in log10 scale
void HistBase::ToGraph(const TString fintxt, const Int_t istart, const Int_t nvar, TGraphAsymmErrors * gr[], const Double_t norix[][4], const Bool_t kpos, const Bool_t kconvert, const Double_t *ymax, const Double_t *ymin, const Bool_t cutyerr, const Bool_t kxerr)
{
  ifstream fr(fintxt);
  if(!fr){
    printf("HistBase::ToGraph file not open!! %s\n", fintxt.Data());
    exit(1);
  }

  const Int_t maxnpnt = 100;
  Double_t xx[maxnpnt];
  Double_t exl[maxnpnt];
  Double_t exh[maxnpnt];
  Double_t yy[nvar][maxnpnt];
  Double_t yerr[nvar][maxnpnt];
  Int_t npnt = 0;

  Char_t buff[500];
  while(fr.getline(buff, 500)){
    TString pnum(buff);
    if(pnum.Contains("#"))
      continue;

    //---------- remove words in front of istart ------------
    for(Int_t ii=0; ii<istart; ii++){
      pnum = pnum( pnum.First(" ")+1, pnum.Length());
    }

    //---------- read x, remove 1 word ----------
    sscanf(pnum.Data(),"%lf",&(xx[npnt]));
    pnum = pnum( pnum.First(" ")+1, pnum.Length());

    if(kxerr)
      sscanf(pnum.Data(), "%lf", &(exl[npnt]));
    else
      exl[npnt]=0;
    pnum = pnum( pnum.First(" ")+1, pnum.Length());

    if(kxerr)
      sscanf(pnum.Data(), "%lf", &(exh[npnt]));
    else 
      exh[npnt] = 0;
    pnum = pnum( pnum.First(" ")+1, pnum.Length());

    if(exl[npnt]<0||exh[npnt]<0){
      printf("HistBase::ToGraph blocking negative errors: %f %f %f\n", xx[npnt], exl[npnt], exh[npnt]);
      continue;
    }

    if(kconvert){
      const Double_t tmpxl = xx[npnt]-exl[npnt];
      const Double_t tmpxh = xx[npnt]+exh[npnt];
      
      const Double_t bg = pow(10, xx[npnt]);
      xx[npnt] = bg;
      
      const Double_t bl = pow(10, tmpxl);
      exl[npnt]=bg-bl;

      const Double_t bh = pow(10, tmpxh);
      exh[npnt]=bh-bg;
    }
    //---------- for each y, read y and yerr, each remove 1 word ----------
    Bool_t kpass=kTRUE;
    for(Int_t ivar=0; ivar<nvar; ivar++){
      sscanf(pnum.Data(), "%lf", &(yy[ivar][npnt]));
      pnum = pnum( pnum.First(" ")+1, pnum.Length());

      //----------- get the absolute values -------------
      if(kpos){
        yy[ivar][npnt] = TMath::Abs(yy[ivar][npnt]);
      }
      if(ymax && ymax[ivar]!=-999 && yy[ivar][npnt]>ymax[ivar]){
        printf("HistBase::ToGraph blocking %s %s -- %d %e > %e\n", fintxt.Data(), buff, ivar, yy[ivar][npnt], ymax[ivar]);
        kpass=kFALSE;
        break;
      }
      if(ymin && ymin[ivar]!=-999 && yy[ivar][npnt]<ymin[ivar]){
        printf("HistBase::ToGraph blocking %s %s -- %d %e < %e\n", fintxt.Data(), buff, ivar, yy[ivar][npnt], ymin[ivar]);
        kpass=kFALSE;
        break;
      }

      sscanf(pnum.Data(), "%lf", &(yerr[ivar][npnt]));
      pnum = pnum( pnum.First(" ")+1, pnum.Length());

      if(cutyerr){
        if(yerr[ivar][npnt]>yy[ivar][npnt]){
          printf("HistBase::ToGraph blocking %s %s -- %d %e < %e\n", fintxt.Data(), buff, ivar, yy[ivar][npnt], yy[ivar][npnt]);
          kpass = kFALSE;
          break;
        }
      }
      
    }
    if(!kpass)
      continue;

    npnt++;
    if(npnt>maxnpnt){
      printf("HistBase::ToGraph npnt>maxnpnt!! %d %d\n", npnt, maxnpnt);
      exit(1);
    }
  }

  //------------- normalize all y to certain y(ipnt), then fill to gr  ------------
  Double_t oldy[nvar][maxnpnt];
  memcpy(oldy, yy, sizeof(yy));

  for(Int_t ivar=0; ivar<nvar; ivar++){
    if(norix && norix[ivar][0]>=0){
      //must use oldy otherwise for ivar > norix[ivar][0]
      Double_t ynorval=-999, ynorerr=-999;
      const Int_t polstart = (const Int_t) (norix[ivar][2]);
      const Int_t poln     = (const Int_t) (norix[ivar][3]);
      MathBase::Polint(poln, &(xx[polstart]), &(oldy[int(norix[ivar][0])][polstart]), norix[ivar][1], ynorval, ynorerr);
      printf("HistBase::ToGraph normalizing y[%d] at %e with %e %e\n", ivar, norix[ivar][1], ynorval, ynorerr);
      for(Int_t ipnt=0; ipnt<npnt; ipnt++){
        yy[ivar][ipnt] /= ynorval;
        yerr[ivar][ipnt] /= ynorval;
      }  
    }

    //--------- Set Maximum of gr -----------
    Double_t maxg = -1e10;
    for(Int_t ipnt=0; ipnt<npnt; ipnt++){
      if(maxg<yy[ivar][ipnt])
        maxg = yy[ivar][ipnt];

      printf("HistBase::ToGraph Please Check %d %d: %e %e, %e %e , %e\n", ivar, ipnt, xx[ipnt], yy[ivar][ipnt], exl[ipnt], exh[ipnt], yerr[ivar][ipnt]);
    }
    //gr[ivar] = new TGraphErrors(npnt, xx, yy[ivar], 0x0, yerr[ivar]);
    //gr[ivar] = new TGraphAsymmErrors(npnt, xx, yy[ivar], 0x0, 0x0, yerr[ivar], yerr[ivar]);
    
    gr[ivar] = new TGraphAsymmErrors(npnt, xx, yy[ivar], exl, exh, yerr[ivar], yerr[ivar]);

    gr[ivar]->SetMaximum(maxg);
    /*
      printf("test....\n");
      for(int ii=0; ii<npnt; ii++)
      printf("%d %d %e %e %e\n", ivar, ii, xx[ii], yy[ivar][ii], yerr[ivar][ii]);
    */
  }
}

TGraph *HistBase::DrawFunc(const TGraph * gr, Double_t (*func)(const Double_t *, const Double_t *), const Double_t *par, const Int_t color,  const Int_t npch0, const Double_t *patch0, const Int_t npch1, const Double_t *patch1)
{
  TString gname(Form("%sfit",gr->GetName()));
  TGraph *gfit = 0x0;
  gDirectory->GetObject(gname, gfit);
  delete gfit;

  gfit = new TGraph(gr->GetN()+npch0+npch1);
  gfit->SetName(gname);

  Int_t ncalc=0;
  for(Int_t ii=0; ii<npch0; ii++){
    const Double_t xx = patch0[ii];    
    const Double_t yy = func(&xx, par);
    gfit->SetPoint(ncalc++, xx, yy);
  }
  for(Int_t ip=0; ip<gr->GetN(); ip++){
    Double_t xx=-999, yy=-999;
    gr->GetPoint(ip,xx,yy);
    yy = func(&xx, par);
    gfit->SetPoint(ncalc++, xx, yy);
  }
  for(Int_t ii=0; ii<npch1; ii++){
    const Double_t xx = patch1[ii];
    const Double_t yy = func(&xx, par);
    gfit->SetPoint(ncalc++, xx, yy);
  }
  printf("npoints:  %d - (%d %d %d) = %d\n", ncalc, npch0, gr->GetN(), npch1, ncalc-npch0-gr->GetN()-npch1);

  gfit->SetMarkerStyle(24);
  gfit->SetMarkerSize(1);
  gfit->SetMarkerColor(color);
  gfit->SetLineColor(color);
  gfit->SetLineWidth(2);

  return gfit;
}

TH1D* HistBase::DrawFunc(const TH1 * hh, Double_t (*func)(const Double_t *, const Double_t *), const Double_t *par, const TString opt, const Int_t color)
{
  TString hname(Form("%sfit",hh->GetName()));
  TH1D *hfit = 0x0;
  gDirectory->GetObject(hname, hfit);
  delete hfit;

  hfit = (TH1D*) hh->Clone( hname );
  for(Int_t ib=1; ib<=hfit->GetNbinsX(); ib++){
    const Double_t ix = hfit->GetBinCenter(ib);
    const Double_t iy = func(&ix, par);
    hfit->SetBinContent(ib, iy);
  }

  hfit->SetLineWidth(2);
  hfit->SetLineColor(color);
  hfit->Draw(opt+" same");

  return hfit;
}

void HistBase::DrawPull(TPad * & ppull, TH1D * & hpull, TF1 * &f0, const TH1 * hh, const TH1 *hfit, const Double_t x0, const Double_t y0, const Double_t x1, const Double_t y1)
{
  //-------------------------------------------------------
  //------------------- Define Pull Hist ------------------

  TString hname(Form("%spull",hh->GetName()));
  gDirectory->GetObject(hname, hpull);
  delete hpull;

  hpull = (TH1D*) hh->Clone(hname);
  for(Int_t ib=1; ib<=hpull->GetNbinsX(); ib++){
    const Double_t cont = hh->GetBinContent(ib);
    const Double_t iy = hfit->GetBinContent(ib);
    const Double_t ip = cont? (iy - cont)/hh->GetBinError(ib) : -999;
    hpull->SetBinContent(ib, ip);
  }

  hpull->SetTitle("");

  //hpull->SetYTitle("Deviation in #sigma_{stat}");
  hpull->SetYTitle("(Fit-Data) /#sigma_{stat}");
  hpull->GetYaxis()->SetTitleSize(0.13);
  hpull->GetYaxis()->SetTitleOffset(0.6);
  hpull->GetYaxis()->SetLabelSize(0.08);
  hpull->GetYaxis()->SetLabelOffset(0.02);
  hpull->GetYaxis()->SetTickLength(0.02);

  hpull->SetXTitle(hh->GetXaxis()->GetTitle());
  hpull->GetXaxis()->SetTitleSize(0.09);
  hpull->GetXaxis()->SetTitleOffset(0.9);
  hpull->GetXaxis()->SetLabelSize(0.07);
  hpull->GetXaxis()->SetLabelOffset(0.02);
  hpull->GetXaxis()->SetTickLength(0.02);
  hpull->SetMaximum(4.5);
  hpull->SetMinimum(-4.5);

  //-------------------------------------------------------
  //------------ Define Pad, draw Pull Hist ---------------

  //TPad(const char* name, const char* title, Double_t  xlow, Double_t  ylow, Double_t  xup, Double_t  yup, Color_t  color = -1, Short_t  bordersize = -1, Short_t  bordermode = -2)
  TString pname(Form("%sppull", hh->GetName()));
  gDirectory->GetObject(pname, ppull);
  delete ppull;

  ppull = new TPad(pname, "", x0, y0, x1, y1);
  //IniDraw(ppull,-999,-999);
  ppull->SetFillStyle(0);
  ppull->Draw();
  ppull->cd();
  ppull->SetBottomMargin(0.17);
  ppull->SetTopMargin(0);
  ppull->SetRightMargin(0.03);
  ppull->SetLeftMargin(0.18);
  ppull->SetFillColor(0);
  ppull->SetTicks(1,1);

  //ppull->SetGrid();
  hpull->SetMarkerStyle(24);
  hpull->SetMarkerSize(0.5);
  hpull->Draw("e");

  //-------------------------------------------------------
  //--------------- Define 0-function, draw ---------------
  TString fname("f0");
  gDirectory->GetObject(fname, f0);
  delete f0;

  f0 = new TF1(fname, "0", -1e10, 1e10);
  f0->SetLineWidth(1);
  f0->SetLineColor(kBlack);
  f0->Draw("same");
}

void HistBase::DrawPar(TPad * &ppar,  TPaveText * &pt, const TObject * hh, const TString sfunc, const Int_t npar, const TString spar[], const Double_t par[], const Double_t err[], const Bool_t kpos, const Double_t px, const Double_t py)
{
  TString pname(Form("%sppar",hh->GetName()));
  gDirectory->GetObject(pname,ppar);
  delete ppar;

  ppar = new TPad(pname, "", 0.699-px, 0.5-py, 0.95-px, 0.91-py);
  //ppar = new TPad(pname, "", 0.599-px, 0.5-py, 1.05-px, 0.91-py);
  //IniDraw(ppar,-999,-999);
  ppar->SetFillStyle(0);
  ppar->Draw();
  ppar->cd();
  
  ppar->Range(0,0,1,1);//redefine pad coordinate system                                                                               
  pt = new TPaveText(0,0,1,1);
  pt->SetBorderSize(-1);
  TText *tt = 0x0;
  const Double_t ts = 0.09;

  pt->SetFillStyle(0);
  tt = pt->AddText(sfunc);
  tt->SetTextSize(ts*0.7/(TMath::Max(sfunc.Length()/74.,1.)));
  if(err[0]>=0){
    tt = pt->AddText(Form("#chi^{2}/NDF %.1f/%.0f", par[npar], par[npar+1]));
    tt->SetTextSize(ts);
  }
  for(Int_t ipar = 0; ipar<npar; ipar++){
    //tt = pt->AddText(Form("p_{%d} %.3e #pm %.1e", ipar, fPar[ipar], fErr[ipar]));                                             

    const Double_t pval = kpos? fabs(par[ipar]) : par[ipar];      
    TString mline;
    if(err[ipar]>0)
      mline = (Form("%s %.3e #pm %.1e", spar[ipar].Data(), pval, err[ipar]));
    else
      mline = (Form("%s %.3e", spar[ipar].Data(), pval));
    tt = pt->AddText(mline);
    tt->SetTextSize(ts);
  }
  pt->Draw();
}

void HistBase::NormalizeSlicesX(TH2D * hh, const Int_t nthres)
{
  Double_t ntot=0;
  Double_t hmin=1e10;
  Double_t hmax=-1e10;
  for(Int_t ix=0; ix<=hh->GetNbinsX()+1; ix++){
    const Double_t xx = hh->GetXaxis()->GetBinCenter(ix);
    TH1D * tmphi = GetSliceX(hh,xx, Form("test%d",ix));
    if(tmphi->GetEntries()<1){
      delete tmphi;
      continue;
    }

    ToPDF(tmphi, 0);
    ntot+=tmphi->GetEntries();

    for(Int_t iy=0; iy<=hh->GetNbinsY()+1; iy++){
      const Double_t cont = (tmphi->GetEntries()<nthres) ? 0 : tmphi->GetBinContent(iy);
      if(cont!=0){
        if(cont>hmax)
          hmax=cont;
        if(cont<hmin)
          hmin=cont;
      }
      const Double_t err = tmphi->GetBinError(iy);
      hh->SetBinContent(ix,iy, cont);
      hh->SetBinError(ix,iy,err);
    }

    delete tmphi;
  }
  hh->SetEntries(ntot);
  hh->SetMinimum(hmin*0.999);
  hh->SetMaximum(hmax*1.001);
}

TGraphAsymmErrors * HistBase::Hist2Graph(const TH1 *hist, const Bool_t delog, const TString xtit)
{
  static const Int_t nmax = hist->GetNbinsX();

  Double_t xx[nmax];
  Double_t errlx[nmax];
  Double_t errhx[nmax];
  Double_t yy[nmax];
  Double_t erry[nmax];

  Int_t ip=0;
  for(Int_t ib=0; ib<=nmax+1; ib++){
    const Double_t cont = hist->GetBinContent(ib);
    if(!cont)
      continue;

    xx[ip] = hist->GetBinCenter(ib);
    errlx[ip] = errhx[ip] = hist->GetBinWidth(ib)/2.;

    if(delog){
      const Double_t x1 = pow(10, xx[ip] - errlx[ip]);
      const Double_t x2 = pow(10, xx[ip] + errhx[ip]);
      xx[ip] = pow(10, xx[ip]);
      errlx[ip] = xx[ip]-x1;
      errhx[ip] = x2-xx[ip];
    }

    yy[ip] = cont;
    erry[ip] = hist->GetBinError(ib);
    ip++;
  }


  TGraphAsymmErrors *gr=new TGraphAsymmErrors(ip, xx, yy, errlx, errhx, erry, erry);
  gr->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
  gr->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
  gr->SetTitle(hist->GetTitle());
  gr->SetMinimum(hist->GetMinimum());
  gr->SetMaximum(hist->GetMaximum());

  const Double_t xmin = hist->GetXaxis()->GetXmin();
  const Double_t xmax = hist->GetXaxis()->GetXmax();
  if(delog){
    gr->GetXaxis()->SetLimits(pow(10,xmin), pow(10,xmax));
  }
  else{
    gr->GetXaxis()->SetLimits(xmin, xmax);
  }

  if(xtit.Length())
    gr->GetXaxis()->SetTitle(xtit);

  return gr;
}

TGraphErrors *HistBase::ProjectRMS(const TH2D * h2, const Bool_t xlog)
{
  const TAxis *ax=h2->GetXaxis();
  const Int_t nx = ax->GetLast()-ax->GetFirst()+1;
  TGraphErrors *gr=new TGraphErrors(nx);
  for(Int_t ix=ax->GetFirst(); ix<=ax->GetLast(); ix++){
    const Double_t xx= h2->GetBinCenter(ix);
    const TH1D * slix = GetSliceX(h2,xx);

    const Double_t px = xlog? pow(10, xx) : xx;
    const Double_t py = slix->GetRMS();
    const Double_t pey = slix->GetRMSError();

    gr->SetPoint(ix-1, px, py);
    gr->SetPointError(ix-1, 0, pey);
    delete slix;
  }

  return gr;
}

#endif
